#######################################################
## 					REPLICATION MATERIAL 
## 
## "Who's to Blame: Political Centralization and 
##  Electoral Punishment under Authoritarianism"
##
##
##					   	   by 				 			 	  
##		Quintin H. Beazer & Ora John Reuter
##
##						     in							 	  
## 		The Journal of Politics		 	  
##													 	  
#######################################################


# File created and checked Dec 2017 using R 3.3.3 (x64). To ask questions 
# or report a suspected error, please contact the authors. We are happy 
# to correspond with others about the project. Send emails to: 
# Quintin Beazer (qbeazer@fsu.edu) or John Reuter (reutero@uwm.edu) 


## TABLE OF CONTENTS
#####################

# This .R script contains code for reproducing the figures  
# presented within the article's main body and in the appendix.
# There is also code for Table A10, the sensitivity analysis results.
# Replication code for the main analyses appears separately in the 
# accompanying .do files.  This file will help reproduce the following:



## FIGURES
# Figure 1: Geographic Variation in the Distribution
#              of Appointed Mayors, 2003-2012

# Figure 2: Marginal Effects of Increasing Unemployment 
#           on United Russia's Vote Share, 
#           Conditional on Mayoral Appointments

# Figure A2: Marginal Effects of Appointment 
#            on United Russia's Vote Share, 
#            Conditional on Changes in Unemployment

# Figure A3: Marginal Effects of Increasing Unemployment Levels
#           on United Russia's Vote Share, 
#           Conditional on Mayoral Appointments

# Figure A4: Balance Tests: Before Reforms, 
#            Non-Appointment Cities & Pre-Appointment 
#            Cities Similar on Most Observable Characteristics


## TABLE
# Table A10: Sensitivity Analysis: Marginal Effects Estimates 
#            of Increasing Unemployment in Appointment Cities, 
#            Depending on Assumptions about Unobserved Confounders


#=============================================
#   PREAMBLE                              ####
#=============================================
graphics.off() # CLOSE ALL OPEN GRAPHICS WINDOWS
rm(list = ls())   # REMOVE ALL OBJECTS


# ----------------- LOAD PACKAGES ---------------- #

## to install packages, use: install.packages("XXXX")
library(foreign)
library(tidyverse)
library(ggplot2)
library(rgdal)



# ----------------- FILE PATHS ---------------- #

## NOTE: Replace the placeholder XXXXX in the file path below 

path <- "C:/XXXX/XXXX/" # PROJECT LOCATION SET BY USER

setwd(path)






##############################################
##   FIGURE 1:  
##   Geographic Variation in the Distribution
#              of Appointed Mayors, 2003-2012 
##   
##############################################


#### load data & create new variables ####


# city data
cdata <- read.dta("RusCityMap.dta")

# shape files of Russia
area <- readOGR(dsn ="RUSmap", layer = "RUS_adm1")
head(area@data)
area@bbox

# turning into dataframe for ggplot (uses broom package)
area <- broom::tidy(area)             # uses command from broom package in tidyverse
area2 <- area %>% filter(long>20)      # to make map more compact


## label for appointed/elected cities
cdata$app <- ifelse(cdata$appcity==1,"Appointed", "Elected")



## set base plot
basePlot  <- ggplot()
basePlot <- basePlot + 
  geom_polygon(data=area2, 
               aes(x = long, y = lat, group=group),
               colour="black", fill="white")


## alter plot

set.seed(84015)
jitPlotbw <- basePlot + 
  geom_jitter( data=cdata, alpha=0.75, 
               position=position_jitter(width=0.5, height=0.5),
               aes(x=lon, y=lat, size = pop, color=app)) +
  theme_classic() + 
  theme(axis.ticks = element_blank(),
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 14, face = "bold"),
        legend.key = element_rect(colour = "white"),
        legend.justification=c(1,0), 
        legend.position=c(1,0)) +
  scale_size(name="Population",
             breaks=c(1e+05, 1e+06),
             labels=c("100K", "1M"))+
  scale_colour_grey(name="", end = 0.5,
                    breaks=c("Appointed", "Elected"),
                    labels=c("Appointed", "Elected")) +
  guides(colour = guide_legend(override.aes = list(size=4)))


## printing
windows(height = 9, width = 18)
jitPlotbw

## saving file in various formats
# png(height = 9, width = 18, 
#     units = "in", res = 800,
#     filename = "AppMap.png")

# tiff(height = 9, width = 18, 
#     units = "in", res = 300,
#     filename = "AppMap.tif")

# jitPlotbw
# dev.off()





## MARGINAL EFFECTS PLOTS
#########################

## load data

mdata <- read.dta("mayorsJOP.dta")



#####################################################################
### FIGURE 2: Marginal Effects of Increasing Unemployment 
### on United Russia's Vote Share, Conditional on Mayoral Appointments 
#####################################################################

# using values from Table 1, Model 3

beta1 <- .0158673 	## appointed
beta2 <- .8840279 ## diff in unemp
beta3 <-  -3.845151

cov13 <- .5805311
cov23 <-  -.31649136
var1 <-  4.6564748   
var2 <- .37834766 
var3 <-  2.4651127

x1 <- mdata$appointed
x2 <- mdata$diff_unemp[mdata$diff_unemp<8 & mdata$diff_unemp>-8]
x2.min <- min(x2, na.rm=T)
x2.med <- median(x2, na.rm=T)
x2.max <- max(x2, na.rm=T)
x1.min <- min(x1, na.rm=T)
x1.med <- median(x1, na.rm=T)
x1.max <- max(x1, na.rm=T)



dx2.x1lo <- beta2 + beta3*x1.min
dx2.x1hi <- beta2 + beta3*x1.max


temp.lo <- beta2 + beta3*x1.min
temp.high <- beta2 + beta3*x1.max

dx2.min.lo <- beta2 + beta3*x1.min - 1.96*sqrt(var2 + (x1.min^2)*var3 + 2*x1.min*cov23)
dx2.max.lo <- beta2 + beta3*x1.max - 1.96*sqrt(var2 + (x1.max^2)*var3 + 2*x1.max*cov23)

dx2.min.high <- beta2 + beta3*x1.min + 1.96*sqrt(var2 + (x1.min^2)*var3 + 2*x1.min*cov23)
dx2.max.high <- beta2 + beta3*x1.max + 1.96*sqrt(var2 + (x1.max^2)*var3 + 2*x1.max*cov23)


## setting ranges for x and y axes

z <- c(0,1)
y <- c( 5, -10 )


## plotting
windows(height=7, width=7)

#pdf("MEunemployment.pdf", height=7, width=7)

#tiff(filename = "MEunemployment.tif", height=7, width=7, 
#    units = "in", res = 300)
    
par(mar=c(4.5,6,2,1))

plot(z,y,type="n", xaxt="n" , bty="n",
     las=1, cex.lab = 1.4, cex.axis = 1.4,
     xlab="",
     ylab= expression("Marginal Effect of" ~Delta~"Unemployment"),
     mgp=c(4,1,0)
)
axis(1, at=c(.3,.7), label=c("Elected","Appointed"), cex.axis=1.4, 
     tick=F,las=1, mgp=c(4.5,1,0))	

segments(0.3,dx2.min.lo,0.3,dx2.min.high, lwd=4, col="blue")
segments(0.7,dx2.max.lo,0.7,dx2.max.high, lwd=4, col="red")
points(c(0.3,0.7), c(mean(dx2.x1lo),mean(dx2.x1hi)), 
       pch=21, col="black", bg="white", cex=2.25)
abline(h=0,lty=2, col="grey70")

#dev.off()




#####################################################################
### FIGURE A2: Marginal Effects of Appointment 
### on United Russia's Vote Share, Conditional on Changes in Unemployment
#####################################################################

# using values from Table 1, Model 1 (in code segment above)

dx1.x2lo <- beta1 + beta3*x2.min
dx1.x2hi <- beta1 + beta3*x2.max

windows(height=7, width=9)
par(mar=c(5,5,2,2))


#pdf("MEappointment.pdf", height=7, width=9)

#png(filename = "MEappointment.png", height=7, width=9, 
#    units = "in", res = 600)

## setting ranges for x and y axes
z <- c(x2.min,x2.max)
y <- c( min(c(dx1.x2lo,dx1.x2hi)), max(c(dx1.x2lo,dx1.x2hi)) )

## plotting
plot(z,y,type="n",xlim=c(x2.min,x2.max),
     ylim= c( min(c(dx1.x2lo-5,dx1.x2hi-5,0)), max(c(dx1.x2lo+5,dx1.x2hi+5)) ), #axes=F,
     las=1, cex.lab = 1.4,
     xlab=expression("Unemployment"[t]-"Unemployment"[t-1]),
     ylab="Marginal Effects of Appointment (dy/dx)"
)

abline(h=0, lty=2, col="grey50")
segments(x2.min,dx1.x2lo,x2.max,dx1.x2hi,lwd=4)
curve(beta1 + beta3*x - 1.96*sqrt(var1 + (x^2)*var3 + 2*x*cov13), lty=1, add=T)
curve(beta1 + beta3*x + 1.96*sqrt(var1 + (x^2)*var3 + 2*x*cov13), lty=1, add=T)

#dev.off()





#####################################################################
### FIGURE A3: Marginal Effects of  Unemployment Levels
### on United Russia's Vote Share, Conditional on Mayoral Appointments 
#####################################################################

# using values from Table 1, Model 3

beta1 <- 8.617566	  ## appointed
beta2 <-  -.0447584 ## unemployment rate
beta3 <-  -5.965547

cov13 <-  -4.1195528
cov23 <- -.77888269
var1 <-   9.2762322   
var2 <-   1.0577606
var3 <-    3.4404521 

x1 <- mdata$appointed
x2 <- mdata$unemp
x2.min <- min(x2, na.rm=T)
x2.med <- median(x2, na.rm=T)
x2.max <- max(x2, na.rm=T)
x1.min <- min(x1, na.rm=T)
x1.med <- median(x1, na.rm=T)
x1.max <- max(x1, na.rm=T)

dx2.x1lo <- beta2 + beta3*x1.min
dx2.x1hi <- beta2 + beta3*x1.max

dx2.min.lo <- beta2 + beta3*x1.min - 1.96*sqrt(var2 + (x1.min^2)*var3 + 2*x1.min*cov23)
dx2.max.lo <- beta2 + beta3*x1.max - 1.96*sqrt(var2 + (x1.max^2)*var3 + 2*x1.max*cov23)

dx2.min.high <- beta2 + beta3*x1.min + 1.96*sqrt(var2 + (x1.min^2)*var3 + 2*x1.min*cov23)
dx2.max.high <- beta2 + beta3*x1.max + 1.96*sqrt(var2 + (x1.max^2)*var3 + 2*x1.max*cov23)


## setting ranges for x and y axes
z <- c(0,1)
y <- c( 5, -10 )


## plotting 
windows(height=7, width=7)

#pdf("MEunemployLevel.pdf", height=7, width=7)

#png(filename = "MEunemployLevel.png", height=7, width=7, 
#    units = "in", res = 600)

par(mar=c(4.5,6,2,1))

plot(z,y,type="n", xaxt="n" , bty="n",
     las=1, cex.lab = 1.4, cex.axis = 1.4,
     xlab="",
     ylab= expression("Marginal Effect of Unemployment"),
     mgp=c(4,1,0)
)
axis(1, at=c(.3,.7), label=c("Elected","Appointed"), cex.axis=1.4, 
     tick=F,las=1, mgp=c(4.5,1,0))	

segments(0.3,dx2.min.lo,0.3,dx2.min.high, lwd=4, col="blue")
segments(0.7,dx2.max.lo,0.7,dx2.max.high, lwd=4, col="red")
points(c(0.3,0.7), c(mean(dx2.x1lo),mean(dx2.x1hi)), 
       pch=21, col="black", bg="white", cex=2.25)
abline(h=0,lty=2, col="grey70")

#dev.off()







#####################################################################
### FIGURE A4: 
##    Balance Tests: Before Reforms, Non-Appointment Cities & 
##    Pre-Appointment Cities Similar on Most Observable Characteristics
#####################################################################


## CREATING THE DATA FRAME FOR PLOT ####

vdata <- matrix(ncol = 2, nrow = 17) %>% 
  data.frame() %>% 
  setNames(c("pval","sig"))

rownames(vdata) <- c("Regional Vote", 
                     "National Vote", 
                     "UR Mayor", 
                     "Opposition Mayor", 
                     "Mayor Age", 
                     "Civil Society", 
                     "Regional Democracy", 
                     "Press Freedom", 
                     "Birth Rate", 
                     "Regional Prominence", 
                     "Caucasus Region", 
                     "Republic Status", 
                     "Russian (%)", 
                     "Unemployment (level)", 
                     "Unemployment (change)", 
                     "Average Salary", 
                     "Working Population")



## T-TESTS ####

## removing obs to compare only elected mayors
bdata <- mdata %>% 
  filter(appointed==0)

## comparing non-appointment to pre-appointment cities
vdata["Regional Vote",] <- t.test(ldvReg ~ appcity, data = bdata)$p.value
vdata["National Vote","pval"] <- t.test(URduma03 ~ appcity, data = bdata)$p.value
vdata["UR Mayor","pval"] <- t.test(URmember ~ appcity, data = bdata)$p.value
vdata["Opposition Mayor","pval"] <- t.test(opp_mayor ~ appcity, data = bdata)$p.value
vdata["Mayor Age","pval"] <- t.test(age ~ appcity, data = bdata)$p.value
vdata["Civil Society","pval"] <- t.test(civsoc91 ~ appcity, data = bdata)$p.value
vdata["Regional Democracy","pval"] <- t.test(reg_democracy ~ appcity, data = bdata)$p.value
vdata["Press Freedom","pval"] <- t.test(as.numeric(press) ~ appcity, data = bdata)$p.value
vdata["Regional Prominence","pval"] <- t.test(pctRegPop ~ appcity, data = bdata)$p.value
vdata["Birth Rate","pval"] <- t.test(birthrate ~ appcity, data = bdata)$p.value
vdata["Caucasus Region","pval"] <- t.test(kavkaz ~ appcity, data = bdata)$p.value
vdata["Republic Status","pval"] <- t.test(republic ~ appcity, data = bdata)$p.value
vdata["Russian (%)","pval"] <- t.test(pctRussian ~ appcity, data = bdata)$p.value
vdata["Unemployment (level)","pval"] <- t.test(unemp ~ appcity, data = bdata)$p.value
vdata["Unemployment (change)","pval"] <- t.test(diff_unemp ~ appcity, data = bdata)$p.value
vdata["Average Salary","pval"] <- t.test(avgSalary ~ appcity, data = bdata)$p.value
vdata["Working Population","pval"] <- t.test(workpop ~ appcity, data = bdata)$p.value


vdata$sig <- ifelse(vdata$pval<=0.05,1,0)
vdata$sig2 <- ifelse(vdata$sig==1,"Signif", "Not")
vdata$varname <- rownames(vdata)




###################
## ## FIGURE ## ##
###################


## window for plotting
windows(width=9,height=6)

## to plot data in ordered manner
vdata2 <- vdata
vdata2$varname <-factor(vdata$varname, levels=vdata[order(vdata$pval), "varname"])

## plotting
dotPlot <- ggplot()

dotPlot <- dotPlot + geom_point(data = vdata2, 
                                aes(x = pval, y = varname, colour=sig2),size=3) +
  scale_colour_manual(values = c("#0277BB","#E65100"), 
                      name="", breaks="", labels="") +
  scale_size(guide=FALSE) +
  xlab("p-value") + ylab("") +
  geom_vline(xintercept=0.05) # + theme_classic() 

dotPlot

# ggsave(
#   "balanceJOP.png",
#   #dotPlot,
#   width = 9,
#   height = 6,
#   units = "in",
#   dpi = 600
# )






########################################################
##  Table A10: Sensitivity Analysis: Marginal Effects Estimates 
#            of Increasing Unemployment in Appointment Cities, 
#            Depending on Assumptions about Unobserved Confounders
########################################################


## based on VanderWeele (2011)


## making a matrix of sensitivity results

  s.lm <- function(the.coef) 
				{
                      gamma=c(.5,1,2,3,5,10)  #effect of latent factor on Y
			    delta=c(.1,.3,.5)
			    bias.mat <- delta %*% t(gamma)
			    colnames(bias.mat) <- gamma
			    row.names(bias.mat) <- delta

                    
                        if(the.coef>=0){
                        ###SENSITIVITY
                        corrected <-the.coef-bias.mat
            }
                        else{
                         ###SENSITIVITY
                        corrected<-the.coef+bias.mat
  }
  return(corrected)
  }

## values Table 1, Model 2
# marginal effects of change in unemployment, given appointments:
# ME = -2.99, CI = [-5.4,-0.57]

s.lm(-2.99)


# (not reported to save space)
# marginal effects of unemployment levels, given appointments:
# ME = -6.18, CI = [-10.24,-2.12]

s.lm(-6.18)



